## Brock, Clare and Daniel J. Mallinson
## Reproduction Code for "Measuring the Stasis: 
## Punctuated Equilibrium Theory and Partisan Polarization"
## 2023-08-23

rm(list=ls()) #clear workspace

library(foreign)
library(moments)
library(lmom)
library(dplyr)
library(tidyr)
library(aTSA)
library(forecast)
library(stats)
library(lmtest)
library(AER)
library(MASS)
library(zoo)

budget <- read.csv("Budget_Authority_Adjusted_2018_for_PAP_1.csv")

budget <- budget[which(is.na(budget$pctChng)==FALSE),]

#budget <- budget[which(budget$majorfunction!=150),] #Uncomment to remove foreign spending

majors <- budget[which(budget$IsItMajorFunction==1),]

pl <- read.csv("US-Legislative-public_laws_20.1_3.csv")

pl <- pl[which(pl$filter_commemorative != 1),]
pl$counter <- 1
pl <- pl[c("id", "year", "pap_majortopic", "counter")]
pl <- aggregate(counter~year + pap_majortopic, data=pl, FUN=sum)

allyears <- data.frame(unique(pl$pap_majortopic), 1948)
names(allyears) <- c("pap_majortopic", "year")

for(i in 2:length(unique(pl$year))){
  year <- unique(pl$year)[i]
  codes <- unique(pl$pap_majortopic)
  data <- data.frame(codes, year)
  names(data) <- c("pap_majortopic", "year")
  allyears <- rbind(allyears, data)
}

pl <- merge(pl, allyears, by=c("year", "pap_majortopic"), all.y=TRUE)
pl$counter[which(is.na(pl$counter)==TRUE)] <- 0

specify_decimal <- function(l,m) format(round(l,m), nsmall=m)

## Kurtosis Change in budget for 5-Year Windows

kurtosis.all <- NULL
majors <- majors[which(is.na(majors$pctChng)==FALSE),]
years <- 1948:2017
for(i in 1:(length(years)-4)){
  use.years <- years[i]:years[i+4]
  use.data <- majors[which(majors$year %in% use.years),]
  kurtosis.all[i] <- round(samlmu(use.data$pctChng)[4], digits=2)
}

## Kurtosis Change in public laws for 5-Year Windows

kurtosis.pl <- NULL
years.pl <- 1948:2020
for(i in 1:(length(years.pl)-4)){
  use.years <- years.pl[i]:years.pl[i+4]
  use.data <- pl[which(pl$year %in% use.years),]
  kurtosis.pl[i] <- round(samlmu(use.data$counter)[4], digits=2)
}

## Polarization Data
#https://voteview.com/articles/party_polarization

pol_dat.all <- read.csv("https://voteview.com/static/articles/party_polarization/voteview_polarization_data.csv")
pol_dat <- pol_dat.all[which(pol_dat.all$year %in% (1948+4):2017),]
pol_dat <- pol_dat[which(pol_dat$chamber=="House"),]

pol.kurt.data <- as.data.frame(cbind(seq(1952,2017,1), kurtosis.all))
names(pol.kurt.data) <- c("year", "kurtosis")
pol.kurt.data <- merge(pol.kurt.data, pol_dat, by="year", all.x=TRUE)
pol.kurt.data$party.mean.diff.d1 <- na.approx(pol.kurt.data$party.mean.diff.d1, method="constant", rule=2)

kurt.pl.data <- as.data.frame(cbind(seq(1952,2020,1), kurtosis.pl))

ts_pol <- ts(data=pol.kurt.data$party.mean.diff.d1, start=1952,end=2017, frequency=1)
ts_kurt <- ts(data=pol.kurt.data$kurtosis, start=1952,end=2017, frequency=1)
ts_kurt.pl <- ts(data=kurt.pl.data$kurtosis.pl, start=1952, end=2020, frequency=1)

png("figure1psj.png", height=5, width=8, units="in", res=600)
plot.new()
plot.window(ylim=c(0,1), xlim=c(1952,2022))
lines(ts_pol, col="red", lwd=1.5)
axis(1, at=seq(1952,2022,10), labels=seq(1952,2022,10))
axis(2, at=seq(0,1,0.1), labels=seq(0,1,0.1), las=2)
#plot.window(ylim=c(0.4,0.8), xlim=c(1952,2020))
lines(ts_kurt, col="blue", lwd=1.5, lty=3)
lines(ts_kurt.pl, col="black", lwd=1.5, lty=2)
title(ylab="Polarization(Red) and Kurtosis (Blue, Black)")
legend(1952,1,legend=c("Polarization", "Budget Kurtosis", "Public Laws Kurtosis"), col=c("red", "blue", "black"), box.lty=0, lty=c(1,3,2))
dev.off()

######### Budget Models- Table 1 ###############

## Kurtosis Change in budget for 10-Year Windows

kurtosis.all <- NULL
majors <- majors[which(is.na(majors$pctChng)==FALSE),]
years <- 1948:2017
for(i in 1:(length(years)-9)){
  use.years <- years[i]:years[i+9]
  use.data <- majors[which(majors$year %in% use.years),]
  kurtosis.all[i] <- round(samlmu(use.data$pctChng)[4], digits=2)
}

## Kurtosis Change in public laws for 10-Year Windows

kurtosis.pl <- NULL
years.pl <- 1948:2020
for(i in 1:(length(years.pl)-9)){
  use.years <- years.pl[i]:years.pl[i+9]
  use.data <- pl[which(pl$year %in% use.years),]
  kurtosis.pl[i] <- round(samlmu(use.data$counter)[4], digits=2)
}

## Polarization Data
pol_dat <- pol_dat.all[which(pol_dat.all$year %in% (1948+9):2017),]
pol_dat <- pol_dat[which(pol_dat$chamber=="House"),]

pol.kurt.data <- as.data.frame(cbind(seq(1957,2017,1), kurtosis.all))
names(pol.kurt.data) <- c("year", "kurtosis")
pol.kurt.data <- merge(pol.kurt.data, pol_dat, by="year", all.x=TRUE)
pol.kurt.data$party.mean.diff.d1 <- na.approx(pol.kurt.data$party.mean.diff.d1, method="constant", rule=2)

kurt.pl.data <- as.data.frame(cbind(seq(1957,2020,1), kurtosis.pl))

ts_pol <- ts(data=pol.kurt.data$party.mean.diff.d1, start=1957,end=2017, frequency=1)
ts_kurt <- ts(data=pol.kurt.data$kurtosis, start=1957,end=2017, frequency=1)
ts_kurt.pl <- ts(data=kurt.pl.data$kurtosis.pl, start=1957, end=2020, frequency=1)

divided.budget <- c(1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,0,0,0,0,1,1,0,0,1,1,1,1,1,1,0)

#Add test for autocorrelation
dwtest(ts_kurt[1:61]~ts_pol)

# LM for Polarization and Budget changes
lm.polar <- lm(ts_kurt[1:61]~ts_pol + divided.budget)
summary(lm.polar)
nobs(lm.polar)

# Test for Stationary
Box.test(ts_kurt[1:61], lag=5, type="Ljung-Box") #non-stationary
Box.test(ts_pol, lag=5, type="Ljung-Box") #non-stationary

# Add time trend
lm.polar.time <- lm(ts_kurt[1:61] ~ ts_pol + divided.budget + c(1:61))
summary(lm.polar.time)
nobs(lm.polar.time)

#LM with lagged Dv
lm.polar.lag <- lm(ts_kurt[1:61] ~ ts_pol + divided.budget + lag(ts_kurt[1:61]))
summary(lm.polar.lag)
nobs(lm.polar.lag)

######### Public Law Models - Table 2 ###############

# Poisson for Public Laws
pl.pos <- pol_dat.all[which(pol_dat.all$year %in% 1948:2020),]
pl.pos <- pl.pos[which(pl.pos$chamber=="House"),]
divided <- c(0,0,0,1,1,1,0,0,0,0,1,1,1,1,0,0,1,1,1,1,1,1,0,1,1,1,1,0,0,1,0,1,1,1,0,1)
pl.pos <- cbind(pl.pos, divided)
mergeyears <- data.frame(1948:2020)
names(mergeyears) <- "year"
pl.pos <- merge(mergeyears, pl.pos, by="year", all.x=TRUE)
pl.pos$party.mean.diff.d1 <- na.approx(pl.pos$party.mean.diff.d1, method="constant", rule=2)
pl.pos$divided[1] <- 1
pl.pos$divided <- na.approx(pl.pos$divided, method="constant", rule=2)
pl.pos <- merge(pl, pl.pos, by=c("year"))
pl.pos <- pl.pos[c("year", "counter", "party.mean.diff.d1", "divided", "pap_majortopic")]
pl.pos$time <- pl.pos$year-1948

pl.ps <- glm(counter~party.mean.diff.d1 + divided, data=pl.pos, family=poisson())
summary(pl.ps)
nobs(pl.ps)
dispersiontest(pl.ps, alternative="greater")

pl.ngb <- glm.nb(counter~party.mean.diff.d1 + divided, data=pl.pos)
summary(pl.ngb)
nobs(pl.ngb)
coeftest(pl.ngb, vcov = vcovCL, cluster=~year) #Standard errors clustered by years
coeftest(pl.ngb, vcov = vcovCL, cluster=~year + pap_majortopic) #Standard errors clustered by year and major topic
est <- cbind(Estimate = coef(coeftest(pl.ngb, vcov = vcovCL, cluster=~year + pap_majortopic)), confint(coeftest(pl.ngb, vcov = vcovCL, cluster=~year + pap_majortopic)))
exp(est)
poldiff <- max(pl.pos$party.mean.diff.d1) - min(pl.pos$party.mean.diff.d1)
poldiff*((1-exp(est)[2,1])*100)

####### Budget Function Changes - Table 3 #############

## Justice Administration
justice <- budget[which(budget$FunctionName=="Justice Administration"),]
round(samlmu(justice$pctChng[justice$year<1990])[4], digits=2)
round(samlmu(justice$pctChng[justice$year>=1990])[4], digits=2)

## Natural Resources
nr <- budget[which(budget$FunctionName=="Natural Resources"),]
round(samlmu(nr$pctChng[nr$year<1990])[4], digits=2)
round(samlmu(nr$pctChng[nr$year>=1990])[4], digits=2)

## Kurtosis by budget category
for(i in 1:(length(unique(majors$FunctionName))-1)){
  print(unique(majors$FunctionName)[i])
  use.data <- majors[which(majors$FunctionName == unique(majors$FunctionName)[i]),]
  print(round(samlmu(use.data$pctChng[use.data$year %in% 1956:1975])[4], digits=2))
  print(round(samlmu(use.data$pctChng[use.data$year %in% 1976:1993])[4], digits=2))
  print(round(samlmu(use.data$pctChng[use.data$year>=1994])[4], digits=2))
}
